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. . . Abstract 

We propose the new exactly solvable pairing model for bosons corresponding to the 
, attractive pairing interaction. Using the electrostatic analogy, the solution of this model 

in thermodynamic limit is found. The transition from the superfluid phase with the 
Bose condensate and the Bogoliubov - type spectrum of excitations in the weak coupling 
| regime to the incompressible phase with the gap in the excitation spectrum in the strong 

C\J ' coupling regime is observed. 

<N 

1. Introduction. 

q | At present time the discrete-state BCS-type [1] pairing models attract much attention 

mainly in connection with the physics of ultra-small metallic grains (for a review see [2]). 
Previously, the discrete-state BCS model was solved by Richardson [3] in the context of nuclear 
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physics. Later the integrability of the model was shown in ref.[4]. The BCS- type exactly 
solvable discrete pairing model for the system of bosons was first considered by Richardson 
[5]. In the continuum limit the condensate fraction and the Bogoliubov type spectrum of the 
low energy excitations (phonons) was obtained. Recently in ref.[6], [7] the pairing model for 
bosons in the context of finite system of bosons confined to a trap was considered. For this 
problem several generalizations of the simplest pairing model analogous to the BCS (fermionic) 
case (for example, see [8]) have been proposed. The remarkable feature of the pairing models 
for the confined bosons [6] is the phenomenon of Bose - condensation. In both cases at 
certain conditions another interesting phenomenon of the condensate fragmentation have been 
observed. From the theoretical point of view, the BCS-type pairing models are of interest due 
to their connection with the generalized Gaudin magnets, Knizhnik-Zamolodchikov equations, 
conformal field theory, and Quantum Inverse Scattering method (for example, see [9]) which 
was studied in a number of papers [10], [11]. From the mathematical point of view the 
models [5], [6] corresponds to the non-compact group SU(1, 1) which is the particular case of 
general non-compact SU(n, m) groups [12]. The construction is equivalent to the non-compact 
SL(2, R) spin chain with different infinite- dimensional representations at different sites. 

In the present paper we propose and solve the simple modification of the model [5] corre- 
sponding to the attractive pairing interaction. Naively, for the attractive pairing interaction 
the ground state energy is a decreasing function of the particle number. However, one can con- 
sider the simple modification of the model [5] which has the correct behaviour of the ground 



1 



state energy both for the finite system and in the continuum limit. Namely, for the infinite - 
volume system of bosons we consider the Hamiltonian 

H = J2 e p n p ~ 9Y,( a t a -p)( a P ,a ~p') + 9'Y,( a t a p)( a t' a p'^ (*) 

p p,p' p,p' 

where g > 0, Bose creation (annihilation) operators corresponding 

to the plane waves with the momentum p a = 2nn a /l (a = x,y,z, I- is the linear size of the 
system) and e p is the dispersion for the free particles (for example, e p = p 2 /2m). For pairing 
models for bosons confined in the external potential the indices in the Hamiltonian (1) should 
represent the states with definite principal quantum number and angular momenta [3]. For 
the system in the thermodynamic limit the sum over momenta p, p' in the second term of 
eq.(l) should be restricted to the values \p\, \p'\ < P, where P is some cutoff P ~ V, and in 
order to have the correct behaviour of the ground state energy in the continuum limit, the 
constant g should be rescaled as g — > g/V, V = / 3 . Although the last term in eq.(l) is nothing 
else but the constant equal to g'N 2 , where iVj, is the fixed number of bosons, the model have 
the correct ground state and in many aspects is a more realistic one in comparison with the 
model with repulsion considered by Richardson [5]. The model (1) can be applied both for the 
finite system of confined bosons and for the system in the thermodynamic limit. Note, that 
if one considers the model (1), as a result of truncation of the initial realistic interaction, in 
general, the terms of both type should be included. Note also that in the real systems like He 4 
the attractive tail of the potential at large distances is always exist. For the finite systems, 
if the pairing interaction is considered as a residual interaction, the coupling constant can be 
of either sign. Previously, the modification of the model (1) for the case of attraction was 
studied numerically for the particular values of the parameters by Dukelsky and Schuck [6]. 

For the model (1) in the limit of the large number of bosons we find the excitation spectrum 
and the occupation probabilities for an arbitrary value of the coupling constant. As a function 
of the coupling constant g we observe the discontinuous transition between the two different 
regimes for the model (1). In the weak coupling regime there is the Bose condensate and 
the Bogoliubov-typc spectrum of excitations. In the strong coupling limit the condensate 
is absent and there is a gap in the excitation spectrum. Qualitatively the results does not 
depend on the spacing and degeneracy of the energy levels and are valid in the limit of the 
large number of particles. The results are compared with the predictions of the mean- field 
theory approach and the Bogoliubov approximation. For the attractive model the naive mean 
- field approximation gives the exact results in thermodynamic limit in the case of the strong 
coupling, while the Bogoliubov approximation is exact in the weak coupling limit in some range 
of density depending on the coupling constant g. We show that the mean-field (variational) 
approach can be modified in order to take into account the Bose condensate and can be used 
for the model (1) to obtain the exact results in the thermodynamic limit in the whole range 
of parameters. 

In Section 2 we review the exact solution of the model, and present the known gener- 
alizations of the model. We show that this class of models for bosons naturally appears in 
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the quasiclassical limit of the algebraic Bethe ansatz transfer matrix. We also present some 
new generalizations of the model [5] which can be useful for studying the superfluidity in the 
framework of this model. In Section 3 we present the solution of the model (1) in thermody- 
namic limit. In Section 4 we compare the exact solution with the predictions of the mean-field 
theory (variational) approach. 



2. Pairing models for bosons. 

The Hamiltonian for the boson pairing model has the form 



L-l 



(2) 



where the coupling constant g is positive for the repulsion. The operators in (2) are defined 
through the pairs of bosonic creation and annihilation operators 0^ at the ith energy 



level 6; 



1, (7 — 1,2). In terms of the pair creation and annihilation operators 

K = <t>U<t>2i, h = <t>U<t>2i, 



the operators B ± are defined as B + = J2ibf, B~ = and n« = riu + n<n [n C i 

For each energy level e; the operators 0^ describe the pair of states differing by the time 
reversal symmetry. For example, for the translationally invariant system the pair creation 
operator has the form 6+ = 0p0t p (zero total momentum). An arbitrary degeneracy Qi of 
each energy level can also be taken into account. For the correct behaviour of the ground 
state one should re-scale g — > g/L, where L is the number of sites, in eq.(2). The rescaled 
value of g will be substituted in the final results throughout the paper. We define the particle 
density p = N^/L, where iVj, is the total number of bosons. The commutational relations for 
the pair creation and annihilation operators have the form 



2''u 



The commutational relations with the operator of the number of bosons rij are: 



± 



These commutational relations are equivalent to the group algebra of the pseudo-spin genera- 
tors for the group SU(1, 1), which differs by the sign for the commutator Sf; S~ from that 
of the SU(2) algebra [12]. As in the case of the conventional discrete BCS-type models, the 
eigenstates can be constructed directly in terms of the operators 



b + 



bt 



We seek for the eigenstates of the Hamiltonian (2) in the form: 

\<f>(t)) = £ + (t 1 )£ + (t 2 )...£ + (t M )W, 



(3) 
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where the state \u) contains only the unpaired states i.e. defined by the conditions: 

bi\v) = {(f>u<t>2i)\v) = 0, 7li\u) = Vi\u), 

where z/j are the (conserved) numbers of the unpaired bosons at the level i. Note that the state 
(3) is degenerate and does not determine the eigenstate completely. One should introduce the 
additional quantum numbers <7j = ±1 such that for each site (nu — n 2i )|z/) = UiV^v). The 
energy does not depend on but the (angular) momentum depends on the quantum numbers 
(Tj. The complete set of states is given by the formula 

b+r°(btr ■ ■ ■ (bt-iT"- 1 ) \("o,«o), (vi,cti), . . . k-i,^-i)>, 

and can be characterized by integer quantum numbers rii, X^n^ = M, instead of the param- 
eters ti (apart from z/«, <jj). Since at g — > the Hamiltonian reduces to the free one each 
eigenstate of (2) at finite g can be characterized by the integers nf\ Yi n i°^ — M, corre- 
sponding to the state at g — 0. In the limit g — > the set of nf^ parameters tj — > We use 
the following commutational relations for the operator S + (t) which can be proved using the 
commutational relations for bf, bi, rii, and different from the formulas for the spin- 1/2 case: 

n,S + (t) = S + (t)n, + -TT^rK, b^ + (t) = S+(t)6, + 7-^(1 + n % ). (4) 

For the Hamiltonian (2) the equations for the parameters tj (3) are obtained in the same way 
as the formulas for the generalized Gaudin magnets (for example, see [5], [10]). The energy 
eigenvalues and the equations for the eigenstates are: 

M Q + v 2 2 

i i=l a l % t a j^i i l j </ 

The total number of bosons equals: N b = J2f=o v i + 2M. Note that for the Hamiltonian 
(2) the number of pairs and the number of the unpaired particles is conserved. Since the 
operator Arii = nu — n,2i commutes with the generators of SU(1, 1) group one can add the 
term J2i hiArii to the Hamiltonian (2) to obtain the model with the external field hi. In 
this case the equations for ti are the same as for the model (2), while the energy equals 
E = J2i( e i + hi(Ji)vi + 2J2iti- The equations (5) are different from the equations for the BCS 
case by the normalization factors and the sign of the second term at the left-hand side. In the 
same way as in ref.[4], the set of commuting operators Hi (i — 1, . . . L) and their eigenvalues 
can be found. In fact, analogously to the case of the SU(2) group, consider the operators: 

^ = -n l + E 7 7 7T (6) 

where we have denoted by (SiSj) = J2 a =x, y ,z ^i^j an d defined 

S+ = ibt, S-=ibi, S' = ^(l + m). (7) 



4 



Note that in terms of initial SU(1, 1) generators bf, bi, 1 + riu + ri2i, the scalar product has 
the form: 

(SiSj) = --(bfb 3 + bp,) + i (l + n t )(l + nj). 

Since the commutational relations for the operators Sf are the same as for the SU(2) group, in 
analogy with the discrete - state BCS- model [4], [6], the operators (6) commute [Hf, H 3 ] = 0. 
At €i = £j the linear combination J2i gives the Hamiltonian (2) while for general q 7^ £j 
we obtain the Hamiltonian depending on the two sets of parameters: 

tf = £e^ + 2£^|W). (8) 

i i<j \^ 

It was noted in ref.[6] that the choice = (e,j) d leads to the realistic model for bosons confined 
in d - dimensional magnetic trap represented by the external harmonic well potential. The 
equations determining the eigenvalues of the Hamiltonian (8) and the operators (6) are given 
by the equations (5) with the parameters replaced by 

Let us comment on the inclusion of the energy level which corresponds to the single Bose 
creation operator 0q~ (p = level in the system with periodic boundary conditions or the 
n = energy level in the spherically symmetric system) i.e. of the non-degenerate energy 
level. One can formally consider the states build up with two Bose creation operators of the 
form (0i"02") n |t / o), where u = 0, 1, and associate with this state the state |z/ + 2n) of u + 2n 
bosons at the energy level 0. The interaction with the other pairs remains the same i.e. of 
the type (<i>t 4>t)(<i>u<foi) (« ^ 0). Thus the energy level n = is considered on equal footing 
with the other energy levels with the exception of the allowed value of the parameter v — 0, 1 
which corresponds to the special type of interaction of pairs with the particles at the energy 
level 0. 

Let us show that the discrete - state bosonic pairing models presented above as well as 
the new models with the interaction of pairs depending on the energy levels, can be obtained 
in the framework of the Quantum Inverse Scattering Method (for example, see [9]). Consider 
the Monodromy matrix defined in the usual way: 

T (t) = K L 

10-^20 • • • Ln , 

where K = diag(e _,?//29 ; e 11 ^ 9 ) is the usual diagonal twist matrix and the Lax operator obeying 
the Yang-Baxter equation is given by 

L i0 (^t) =£ l -t + r ] (aS l ), (9) 

where the operators 5? (a = x,y,z) were defined through the SU (1,1) generators in the 
previous section, a a are the Pauli sigma matrices and ^ are the inhomogeneity parameters. 
Considering the quasiclassical limit rj — > of the transfer matrix Z(t) = Tr (T (t)), we obtain 
the family of the Hamiltonians depending on the spectral parameter t, which commute at 
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different values of the parameters: 

" <<> = h ? (^6) (1 + + 2 § (^W Si ^ (10) 

[if (t); if (t')] = 0, where the expression for the scalar product (SiSj) was presented in the 
previous section. From eq.(10) taking the limits t — > £j or t — > oo the different pairing models 
for bosons can be obtained. The limit i — > & corresponds to the operator ii; (6). The 
eigenvalues of H(t) can be obtained from the eigenvalues of the transfer matrix Z(t). We 
define the monodromy matrix in the following way: 

Tlt\-( m B{t) \ 

and seek for an eigenstates in the form 

\(f>(t)) = B(t 1 )B(t 2 )...B(t M )\v), 

where the reference state \u) is defined by the conditions S~\v) = and n^v) = Vi\v). As in 
the usual SU (2) case we observe that 

C(t)\v) = 0, 

and the state \u) is an eigenvector of A(t) and D(t). The eigenvalues of the operators A(t) 
and D(t) are: 

mw) = n - * + < i + ^)/2) D(t)\u) = n (e« - * - ^(i + ^)/2) k). 

Following the well known procedure we obtain the Bethe ansatz equations: 
eV/g tt ^-q + ^(l + ^)/2 \ = " a-t Q -7y \ 

The corresponding eigenvalue of the transfer matrix Z(t) equals 

m = n ( ^rr ) n(& - * + v(i + ^)/2) + n ( t ~! t +r? ) n& - * - + *«)/2), 



a 



where the two terms corresponds to the operators A(t) and D(t) respectively. Considering 
the quasiclassical limit of the Bethe ansatz equations, one reproduce the equations (5) for the 
eigenstates of the pairing Hamiltonian (2). The eigenvalues of the operators Hi and H (2) can 
be obtained from the last expression for A(t). According to the general procedure [13] one can 
obtain the fundamental R - matrix corresponding to the direct product of two representations 
of the group SU (1,1) (the so called fundamental Lax operator) which will lead to the new 
transfer matrix Z^\t) with the trace over the infinite- dimensional space, commuting with 
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the transfer matrix Z(t) and the Hamiltonians of the new type, which is beyond the scope 
of the present paper. In order to obtain the trigonometric transfer matrix, one should have 
the special quantum group commutational relations [S + ; S~] = sm(2r]S z )/sm(r]), which are 
not fulfilled for the SU(l, 1) generators. However, since the commutational relations for S ± , 
S z coincide with the usual SU(2) algebra, the Gaudin - type Hamiltonians [8], which are 
quadratic in the operators S a , can be obtained in the trigonometric case. Thus all the known 
results, for the Gaudin- type Hamiltonians for the trigonometric and the elliptic cases, can be 
generalized to the case of SU (1,1)- generators, constructed with the help of bosonic operators. 

Let us briefly comment on the possible generalizations of these models. One can use the 
representation of spin-s operators through the single Bose creation and annihilation operators 
0+, 0, [0; 0+] = 1 as S + = 0+(2s - 0+0) 1/2 , S~ = (2s - 0+0) 1/2 0, S z = 0+0 - s, to include 
this spin in the quasiclassical Hamiltonian. In the fermionic case this leads to the generalized 
Dicke model if the limit s — > oo is taken. If one considers the limit £i — > oo for this single site 
in the hyperbolic version of the model (8) one can eliminate the terms, which do not conserve 
the number of bosons and obtain the model describing the interaction of single oscillator with 
the bosonic degrees of freedom: 

H = Ld(f) + (f) + A(0 + 0) I £i n i + CiCTiVi J + Y hiUi + Y ti^i + gB + B~ . 

This Hamiltonian contains a number of free parameters which can be chosen in order to 
get the realistic model. In the sector with the oscillator excited to n-th energy level the 
model is reduced to the boson pairing model with the renormalized coupling constant g. At 
the same time the excitation energy (level spacing) for the oscillator depends on the average 
occupation numbers rii for bosons. In contrast to the same model for fermions, the occupation 
probabilities (rij) can be a small numbers, which allows for the small renormalization of the 
oscillator frequency. This model can be useful for studying the superfluidity without any 
assumptions. 

3. Continuum limit. 

Let us solve the model (2) in the continuum limit. Assuming that the distribution of roots 
ti can be approximated by the continuous density R(t), which is valid for the large number of 
pairs M, we get the following equation: 

fdt'^l = m , f(t) = l-Wj^-i (12) 

Ja I — I g Z a I — € a 

where the integral in a sense of principal value over the support of the function R(t) is implied 
and C a = Q a + u a . According to [5] for the case of repulsion the ground state corresponds 
to the roots ti located at the interval (eo, ei). One can argue that since at g — > the ground 
state corresponds to all ti — > eo and the roots ti cannot cross the values e a for varying g, all 
ti G (eo, ei). Thus one should solve the equation (12) assuming that the support of the function 
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R(t) is the interval (eo, ei). The structure of the ground state for the repulsion and the general 
behaviour [5] of solutions of the equations (5) can be easily seen from the electrostatic analogy. 
Electrostatic analogy for the equations of the type (5) was previously used for the solution of 
the equations for the case of the BCS problem (for example, see [14]). One can consider the 
functional of the roots ti as an energy of charges at the two-dimensional complex plane: 

®(U) = ~J2 C Mti - e«| - 2^]n\U -tj\ + (2/g)J2MU)- 

i,a i<j i 

This energy functional corresponds to the repulsion of unit charges ti and the repulsion of 
the charges with the charges of the same sign C a > placed at the fixed points e Q . The 
external electric field of the magnitude 1/g is applied. The condition of stationary point (not 
minimum) of the functional $(tj) with respect to the positions of the charges U leads to the 
equations (5). It is convenient to imagine the charges e a placed at the y axis of (x, y) plane as 
an energy levels. Then for the case of repulsion the external electric field is directed down, and 
each root gives the contribution to the energy equal to its height. One can see that all roots ti 
are real, since due to the repulsion and the external electric field all the other configurations 
are unstable. Physically the picture is as follows. For each charge the repulsion due to the 
external charges e a below this charge, and the other roots below this charge, produce the 
force directed up. This force is compensated by the other charges above this charge and the 
external electric field directed down. For the ground state the roots U should be placed as 
low as possible. This picture allows one to use the physical intuition to find the solutions for 
the ground and the excited states of the model (2). For instance, the general behaviour of the 
solutions [5] is obvious. 

Here we consider the pairing model (2) for the case of attraction g < or equivalently 
the model (1) for g > 0. It was already mentioned that due to the additional term (1) the 
behaviour of the ground state energy as a function of particle number is correct. In many 
aspects the model (1) is more realistic in comparison with the model with repulsion [5]. For 
example, it has the Bogoliubov-type spectrum of excitations and the Bose condensate which 
varies continuously with the coupling constant from zero at some critical coupling g c to N b at 
g = 0. Later on we will omit the last term in eq.(l) in all the formulas. In the framework of 
electrostatic analogy the case of attraction corresponds to the external electric field directed 
up. Thus for any coupling constant g for the ground state all roots of the equations (5) located 
below the lowest energy level eo = 0. For small \g\ they are close to eo, while for large \g\ they 
are far below the level eo- At the sufficiently small \g\ the density of roots R(t) grows at t — > 
and bounded from below at some fixed point —b (b > 0). 

The simple method to find the solution for R(t) (12) is, using the electrostatic analogy, to 
consider the electric field at the complex plane z produced by the unit charges ti located at 
the interval T = (a, b) of the real axis, the charges e a , and the external electric field: 

h(z) = fdt^- t -f(z) 
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where the discontinuity Ah(t) at T is given by the density of charges R(t): Ah(t) = h(t + 
iO) —h(t — iO) = 2mR(t). The equation (12) takes the form h(t) = 0, where h(t) is an average 
value of the field at both sides of T, and can be represented in the form: 

J_ / dz h^l = f (t) (13) 

2mJcz-t Jy ' y ' 

for t e T, where the contour C encloses the interval (a, b). For the sufficiently small coupling 
constant we use the following ansatz for the electric field h(z) in the complex plane which has 
the branch cut along the interval (a, b) (in this case we take a = and the interval (—6,0), 
b > and use the coupling constant for the attraction g > 0): 




= [L^-Tj)- (14) 

where the function 0(£) can be fixed from the condition for the residues of h(z) at the points 
e a which are equal to —C a /2, 

0(0 = ' P(0, P(0 = ~\ E CJ{i - e a ). 

The constant term in the parenthesis is fixed from the behaviour of the field at infinity, and 
the value of b is determined from the condition / dtR(t) = M. The number of pairs and the 
energy AE = J2i can be represented as an integrals in the complex plane over the contour 
C enclosing the interval T: 

M= - A£ =-£S 2 ^' < 15) 

The integrals can be evaluated by means of deformation of the contour C into the small 
contours around the points e a and the large circle at the infinity The equivalent way to 
find the energy is to substitute the ansatz for h(z) into the equation (13), which after the 
deformation of the contour C allows one to find the function 0(£)> presented above and the 
term 1/g in eq.(14). Using the same formulas for M and E (15), we obtain the following 
equation for the particle number: 

b ( f \ 1/2 

- = N b + L-Y,C a S(e a ), S(0=f^r) , (16) 

which determines the value of the parameter b. In contrast to the case of repulsion apriory we 
did not have any condition, which determines the upper bound for |6|: the support of R(t) is 
not bounded from below for g — > oo. The value of b found from the last equation should be 
substituted to the equation for the energy (15) which takes the form: 

E = ~ E ^ + E CaS(e a )(e a + b/2) - h —. 

a a ^9 
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Using the equation (16) one can represent the last equation in a more convenient form: 



£ = "E e « + E C a E(e a ) - b -(N b + L)- E(e) = ^e(e + b). (17) 

In order to find the excitation spectrum and the occupation probabilities one should calcu- 
late the variation of the energy (17) over the quantum numbers v a and the energy levels e a 
respectively, taking into account the variation of the parameter b eq.(16). Let us note that 
the units for q can be chosen in an arbitrary way, see eq.(5). The possible choice is ei = 1/L, 
such that Lei = 1- In thermodynamic limit there are two parameters - the density p and 
the coupling constant g. For example, one can imagine a one-dimensional lattice model with 
linear dispersion relation and L lattice sites. We will assume the units Lei — 1 and, as an 
example, consider the equal-spacing L level model with Q a = 1 and use the rescaled coupling 
constant g — > g/L va. the final results. We obtain from the equation (16) at v a = the 
following equation for the parameter b: 

b = g(f(b)+p), (18) 

where f(b) is a smooth function which varies from zero to unity for b G (0, oo). For example, 
for the equal-spacing model with Lei = 1 we have 

f{b) = Mn ((1 + VT+b)/Vb) + 1 - VT+b. 

Equation (18) gives the value of b as a function of the parameters g, p. First, consider the 
limit of the small coupling constant j<1, jj5<l, such that g <C gp. According to the last 
formula this limit corresponds to b = gp, and the high density limit p ~ 1/g. In this limit one 
can neglect the last sum in eq.(16) and obtain the excitation spectrum and the occupation 
numbers. Considering the energy (17) as a function of the quantum numbers v a and taking 
into account the variation of the parameter b according to eq.(16), we find the spectrum of 
phonons: 

E(e) = ^e(e + gp). 

As in ref.[5] one can show that the states corresponding to the excitation of pairs have the same 
energy, so that the (bosonic) quasiparticle interpretation of the excited states is true. This 
formula, corresponding to the particular limit g <C gp <C 1, is in agreement with predictions 
of the Bogoliubov approximation. However, in contrast to the repulsion, this spectrum is 
exact for an arbitrary value of the parameter gp with respect to the energy level spacing ei, 
provided the condition g 1 is satisfied. 

Variation of the energy (31) with respect to the parameters e a gives the occupation proba- 
bilities (n a ) which are different for (n ) (condensate) and (rij), i ^ 0, which can be easily seen 
from the electrostatic analogy. In fact, if the parameter gp is not too large, shifting the level 
eo = down will obviously shift the distribution of roots and the lower boundary —b down as 
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a whole, which means the existence of the condensate. Considering the variation bEjbti for 
i 7^ 0, we obtain: 

ei + gp/2 . 

("<> = 7 . - 1, ^0. (19) 



The condensate fraction iV can be evaluated using the equation N = N b — J2i^o( n i)- At 
e i ^ <?P the sum in (19) can be replaced by the integral, which gives: 



N' = N b -N = L(yJl + gp-l), 



The parameter which governs the condensate fraction is g: in the limit considered, N' = 
(gp)L = gN b « N b . 

Next, consider the case b ~ 1. According to eq.(18) it is possible in the two cases: (i) 
g <^ 1 and p ~ 1/g 1; (m) # ~ 1, p ~ 1. In both cases calculating the excitation spectrum 
and the occupation numbers from the equations (16), (17), i.e. taking the variation of the 
energy (17) with respect to v a and e a (taking into account the variation of the parameter b 
according to eq.(16)) we obtain the expressions 

E(e) = ^e(e + b), (n t ) = J ± ^ - 1, i ? 0. 

y/€i{€i + b) 

and the expression for the condensate 

N' = N b -N = L(VT+b- 1). (20) 

Since gp ~ 1, in the case of large density p ^> 1 (case (i)) we will always have N' N b . 
However, in the case (ii), p ~ 1, for the coupling constant g larger than some critical value g c , 
the last equation gives N' > N b . That means that for the sufficiently large coupling constant 
the solution (14) is not correct. 

Below we will show that at g > g c the solution should be modified. We will also show that 
the critical value g c is determined exactly by the condition N'(b) = N b , where b is the solution 
of the equation (18) (i.e. we will show that this condition coincides with the condition (26), 
see below). Here let us present the physical arguments, which show that at large g the new 
phase with the gap in the excitation spectrum should exist. As a limiting case, consider L 
(the large number) levels glued together. In this case the repulsion directed down is strong in 
comparison with the external field directed up and there cannot be roots tj in the vicinity of 
e = 0. Thus the support of R(t) should be located far below e = 0, at the distance of order 
~ gL. The ground state energy will be of order ~ —gLN b /2, and the gap in the excitation 
spectrum will exist. This picture is in agreement with the energy of the one- level model [5] 
with Q replaced by L. So, as a first step, one could solve the one-level model with Q replaced 
by L and v = at g — > oo, which would be the particular case of the general solution. In 
other words, at large g (weak external field) the distribution of charges will be unstable if the 
length |6| ~ gp is much larger than the length Le±. 
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Thus, we seek for solution of the equation (12) with the density support at the interval 
(b, a), a, b < 0. In general one expects that since there is no external charges e a in the vicinity 
of the interval (b, a), the support of R(t), it should be equal to zero at the endpoints. The 
numerical calculations suggest that for \a\ > the function R(t) is, in fact, equal to zero at 
the points a, b. It might seem that the ansatz for h(z) should be chosen in such a way that 
as limiting case a = it would contain the solution for the interval (6, 0) i.e. in the form (14) 
with a 7^ 0. However, we will show that correct solution reproduce eq.(14) at a = 0. One 
can perform the calculations with the field of the type (14) and find that the parameters a, 
b are not completely fixed from the solution itself and one finds a number of functions R(t) 
with different energy, which is not correct, as can be seen from the electrostatic analogy. Thus 
let us find the solution of (12) with the density R(t) equal to zero at the endpoints of the 
interval (6, a), a, b < 0. Since the field h(z) should be a constant at the infinity, we consider 
the following function: 

h(z) = y/(z-b)(z-a) ( [ d£^j , (21) 

where eo = and a, b < 0. The points a, b should be determined from the solution itself. Note 
that there are no poles other than the poles corresponding to the charges e$ in h(z). After 
changing the signs of the parameters a, b, from the equation (13) we find 

1 



Z a 

and simultaneously the condition for the behaviour of the field at the infinity: 

/«0 = -^E^) = --, 

or, equivalently, 

E = =-■ (22) 

a yj(e a + a) (e Q + b) 9 



The first of the equations (15) gives 

a + b 
9 



N b + L-J2C a e a S(e a ), (23) 



where the relation M = (iVj, — u)/2 was used. Substituting the ansatz (21) to the second of 
the equations (15) and using the equation (22) we obtain the energy: 

£ = -E e « + E C a e 2 a S(e a ) + W C a e a S{e a ){a + b) - ^31. 

Taking into account the equation (23) after some algebra this expression can be represented 
in the following form: 

£ = "E e « + E C a E(ea) -(N b + L)^ + ±(b- a) 2 . (24) 

a a 1 ^9 
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Thus the parameters a, b determined from the equations (22), (23) should be substituted 
to the energy (24). If the parameters a, b are fixed, if a ^ 0, the gap in the spectrum of 
excitations will appear and the Bose condensate will be absent. One can further rewrite the 
equations (23), (24) in order to see the similarity with the mean-field (variational) equations 
presented below. Introducing the notations 



p 



the equations (23), (24) for p, A at v a — take the form: 



^(e l + p) 2 -A 2 ' i ^/(ei + /j,) 2 - A 2 9 

E(p, A) = - £ e t + £ ^ + /i) 2 - A 2 - (iV 6 + LV + — . (25) 

The gap in the excitation spectrum equals \/ p 1 — A 2 . The parameters /i, A found from the 
first two of the equations (25) should be substituted into the energy E(p, A) (25). The first 
two of the equations (25) are equivalent to the condition of minimum of the energy E(a, b) 
over the variables a, 6, which in terms of new variables reads 8E/8p = 0, 8E/8A = 0. Thus 
the equations (25) are equivalent to the equations obtained from the mean-field theory (see 
below). The difference of the exact solution with the mean-field approach can appear only in 
the weak coupling regime in presence of the Bose condensate. 

If the minimum of the energy exist, the solution of the equations (25) can be easily found. 
For example, for the equal-spacing L-level model with Le\ = 1, taking the variations of (25) 
over n and A we get the equations presented in the next section. The condition of the existence 
of the solution is 



VV 2 - A 2 = — — - (2C - p(p + 2)) > 0, C 



2(p + 2) v " ' e 2 /9-l 

Separately the parameters //, A can be found from the equations 

(/i 2 - A 2 ) 1/2 = (C 2 - A 2 )/2C, n = (C 2 + A 2 )/2C. 

For a given density p the last equation gives the critical value of the coupling constant g c : 

M = RTT27rt ' (26) 

For g > g c {p) the solution of the equations (25) exist, \a\ > 0, and the gap in the energy 



spectrum V ab = p 2 — A 2 > 0. The Bose condensate is absent in this phase. For g = g c (p) 
we have a = and for g < g c (p) the gap closes and the solution (14) with the Bose condensate 
described above is valid. In fact, let us show that the critical value (26) coincides with the 
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value of g determined by the condition Nb — Nq < Nb in the framework of the solution (14) by 
the equation (18). One observes that eq.(16), (18) can be represented in the following form: 



The cancellation of the first two terms at the right-hand side of this equation is equivalent to 
the condition N' = Nb in the framework of the solution (21), while the last sum equals b/g in 
the framework of the solution (21) at a = 0. In fact, from eq.(25) at a = (/x = A) we obtain 
exactly 2/g = J2i(^i(^i + b))^ 1 / 2 . Therefore, two estimates of the critical point g c found from 
two solutions in the weak and the strong coupling limits are coincide. 

Let us show that at the value a = the density R(t) given by eq.(21), which was equal 
to zero at this point, is reduced to the density in the weak coupling regime (14) which is 
unbounded at t — 0. Due to the term ~ \ jt in the parenthesis of eq.(21) one can rewrite the 
density (21), 

S(e r 



R{t) = y(t-a)(t-b) 



t-€ r 



in the following form: 



1 It- a ( 1 1 ^ S(e a )(e a - a)\ 

if the condition / d£(f)(g) = —1/g (22) is taken into account. The last expression coincides 
with the result obtained from the ansatz of the type (14) if the condition of minimum of the 
energy (22) as a function of a, 6, ^2 i (l/E(e i )) = 2/g is satisfied. However, let us stress that 
the transition between the two phases at the critical point g c {p) is discontinuous. 

Thus, we have shown that the transition from the strong coupling incompressible phase 
with the gap to the phase with the Bose condensate and the Bogoliubov- type spectrum of 
phonons takes place at the coupling constant g = g c (p) (26). At this point the condensate is 
equal to zero, N Q = 0, but at g < g c {p) the condensate increases according to the equation 
(20) until the value N = N b is reached at g — 0. Let us note that qualitatively these results 
are valid for the model with an arbitrary degeneracy of energy levels fl a and an arbitrary level 
spacing. Numerically the dependence g c (p) will have the different form. The limiting case of 
the solution at |^| — > oo coincides with the solution of the one-level model in this limit. 

4. Mean-field solution. 

Here we consider the mean field or variational solution of the pairing model (2) for the 
case of attraction: ^ 

H=Y,e i n i -gB + B-, g > 0. (27) 

i=0 
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Let us describe the mean-field approach for the model (27) and find the range of the parameters 
for which the solution is exact in the thermodynamic limit. The mean-field Hamiltonian has 
the form 

A 2 

Hmf = E( e * + A*H + A + h) - pN b + — , (28) 

i i 9 

where p is the chemical potential and the variational parameter A is real. The expression (28) 
can be considered in a sense of the Hubbard - Stratanovich transformation in the functional 
integral which can also be used to establish the validity of the mean-field theory. For the 
Hamiltonian (27) the mean field theory (28) is equivalent to the variational procedure with the 
trial wave function analogous to the BCS wave function. It is well known that the variational 
solution for the BCS Hamiltonian is exact in the thermodynamic limit (for example, see [15]). 
In contrast to the BCS case, for the bosonic model one has to introduce the chemical potential 
in order to fix the particle number. We show that in some range of parameters, at g > g c (p), the 
variational solution for the bosonic pairing model coincides with the exact solution presented 
above. At g < g c {p) the naive mean field solution is not correct. However, for our model one 
can modify the mean field (variational) approach taking into account the Bose condensation 
to obtain the exact results presented above in the whole range of parameters (except the 
extremely small coupling constant gp ~ ei). 

Each of the quadratic Hamiltonians Hi in the sum (28) can be diagonalized by means of the 
Bogoliubov transformation. For each site i introduce the new Bose creation and annihilation 
operators Xi,2, X12 according to 

4>t = c xt + SX2, <f>t = c xt + s Xi, 
where the coefficients q, are assumed to be real, 

c--s 2 = l, Q = ch(0i), Si = sh(0i). 

The expectation values of the particle number rii and the energy Hi in the ground state are 

(m) = 2s 2 , (Hi) = + ^)(c 2 + s 2 - 1) + A2c lSi . 
The condition of cancellation of the terms X1X2 and xtxt takes the form: 

2 ^ 2 = th 20. = 

Ci + Sj 6i + p 

Thus we obtain the expressions for the energy and the number of particles as a functions of 
the parameter A and the chemical potential p: 

A 2 
9 



E MF (A) = £ ( v /(6, + / i) 2 -A 2 (l + nf) - (e, + /.)) - pN b + — , 



^'wb' 1 ^- 1 '- (29) 
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where the operator nf equals 

n} = XuXu + xtiX2i, Vi<Ji = xtiXu ~ xtiX2i- 

The parameters A and \x should be determined from the condition of minimum of Emf(A) 
(29) with the condition of fixed number of particles N b . From eq.(29) the excitation energy 
Ei = \J (e$ + jj) 2 — A 2 . The ground state corresponds to the quantum numbers nf = 0, or, 
equivalently, to the state |0) x annihilated by the operators Xi,2'- 

Xu\0) x = 0, X2,|0) x = 0, t = 0,...L-l. 

In terms of the initial operators (frf, 0^ this state can be represented as 

|0> x = ne^^|0>, a l = ^ = th( ( p i ), (30) 

i °i 

where |0) is the vacuum with respect to the initial operators: (f>u\0) = 0, 02i|O) = 0. In fact, 
for each i for the state \a) = exp(a0i~02~)|O) one finds 

(0i -a<j)+)\a) =0, \ a ) = e a{ ^rt ) \0). 

Substituting the expressions for the operators <f>i, <p2, one can see that the state (30) is an- 
nihilated by the operators Xi 5 X2 provided the condition a = s/c = th(0) is satisfied. The 
excited states can be constructed from the state (30) by action of the operators xt, xt : 

(xl) ni (xl) n2 \0)x, v=W-n 2 \. 

Let us show that for the model (27) the mean-field theory approach is equivalent to the 
variational procedure with the trial variational wave function of the form (30). Although this 
wave function does not correspond to a definite particle number, it can be fixed in an average 
as in the usual BCS theory, which is justified in the continuum limit. Expectation value of 
the Hamiltonian (27) over the state (30) as a function of the variational parameters (pi takes 
the form: 

E = J2^ + ^(2s 2 l )+g(^2c l s}j -fiN b . (31) 
Taking the variation of (31) with respect to 4>i one finds the equations presented above with 

A =9^2^). 

i 

The average particle number N b = J2i^i- The existence of the condensate means O — ► 00. 
Substituting this value to the right -hand side of eq.(31) and assuming iVo = N b , one finds fi = 
gp/2, and the spectrum E(e) = \Je(e + go) in agreement with the Bogoliubov approximation. 
Thus although the Bogoliubov approximation corresponds to the variational estimate of the 
energy, in general, it does not correspond to the minimum of the energy on the class of the 
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wave functions (30). However we show that in the weak coupling limit the naive variational 
approach does not lead to the correct results while the Bogoliubov approximation gives the 
exact results in the weak coupling limit if the density is not too small (1 <C p <C 1/fiO- 

In the strong coupling regime g > g c (p) the equations (29) coincide with the exact equations 
obtained in section 6. The equation 5E MF (A)/8A = together with the second of the 
equations (29) allows one to find the parameters p, A, the occupation numbers and the gap 
in the energy spectrum. In particular, for the equal-spacing L level model with Le\ = 1 the 
equations take the form: 

In (l + p + ^(1 + p) 2 - A 2 ) - In (/i + yV 2 - A 2 ) = - g , (32) 

p+l = yV+l) 2 -A 2 - vV-A 2 , (33) 

which have the solution found in section 6. The results are in agreement with the exact 
solution. 

Let us see if the exact solution in the weak coupling case g < g c can be obtained in the 
framework of the mean-field (variational) approach. To get the expectation value (no) = N of 
order N b) one should take A = p+5, where S is the small parameter of order 1/L. Substituting 
the value p = A into one of the equations (32), (33) one finds that the result for A contradicts 
the exact solution. Thus the naive mean field approach fails for the region of the parameters 
where the solution of the equations 5Emf/$P = 0, 5Emf/5A = does not exist. To get the 
correct results, one should use the following method. First, substitute the parameter p = A 
into E MF (p,A) — > E MF (A, A) (29). Then the solution of the equation 5E MF /5A = gives 
the results in agreement with exact solution of section 6. In fact, one can see that the variation 
of this function leads to the equation (16), which in the framework of the exact solution was 
used to determine the parameter b. The validity of this method can be shown in the same way 
as for the usual BCS model. In the framework of the functional integral approach the factor 
L (the volume) appears in the exponent in front of the action if the condensate is absent. To 
take into account the condensate one can introduce the 5- function of the form 

5(n (p,A)-N b + N'(A)), 

where the function no(p, A) = (no) is given by eq.(29) for i = and the function N'(A) is 
determined by the sum J2ijto( n i), with (n^) given by eq.(29) with p = A. This factor will give 
p = A with the accuracy of order l/N b and remove the integration over p. If the saddle point 
for the remaining integration over A exist and gives the value N'(A) < N b , which indeed takes 
place, the solution is exact in the thermodynamic limit. The particle number is correctly fixed 
within this approach. The same can also be shown using the trial variational wave function of 
the form \N , 1; . . . <Pl-i), where N and (pi, i ^ are the variational parameters. Thus the 
modified mean-field approach is valid in the whole range of the parameters with the exception 
of the extremely small coupling constant gp ~ ei, when the value of N', the number of particles 
out of the condensate, becomes of order of unity. 
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Conclusion. 

In the present paper we have shown that the discrete-state BCS-type pairing models for 
bosons can be considered as a quasiclassical limit of the eigenvalue problem of the general 
transfer matrix in the framework of the algebraic Bethe ansatz method. We introduced the 
new pairing model for bosons corresponding to the attractive pairing interaction. It was 
shown that the weak coupling phase, g < g c , is characterized by the Bose condensation and 
the Bogoliubov-type spectrum of phonons. In the strong coupling phase at g > g c the Bose 
condensate is absent and there is a gap in the excitation spectrum. Note that the transition 
of this type from the incompressible Mott insulating phase to the superfluid phase is usually 
expected in the Bose Hubbard model. We have shown that naive variational approach is not 
applicable in the weak coupling limit at g < g c , when the condensate fraction exist. However, 
for our model one can modify the variational procedure taking into account the condensate 
fraction to obtain the exact solution in the whole range of parameters. The Bogoliubov 
approximation gives the correct results in agreement with the exact solution in the limit 
g < 1 and 1 <C p <C 1/g, such that the parameter gp <C 1, i.e. when the parameter b = gp 
(see eq.(18)). The proposed model with an attractive pairing interaction can be interesting 
both in the context of applications to the finite systems of the confined bosons and for studying 
the phenomenon of superfluidity in the exactly- solvable model. 
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